directory = "/Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/"
ba_directory = "ba/"
stereo_directory = "stereo/"

!asp_plot \
  --directory $directory \
  --bundle_adjust_directory $ba_directory \
  --stereo_directory $stereo_directory \
  --map_crs EPSG:32604 \
  --dem_gsd 1 \
  --subset_km 1 \
  --add_basemap True \
  --plot_icesat True \
  --plot_geometry True

Processing ASP files in /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/


Reference DEM: /Users/ben/Dropbox/UW_Shean/COP/COP30_utqiagvik_lzw-adj_proj.tif


ASP DEM: /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/stereo/20220417_2252_1040010074793300_1040010075633C00-DEM_1m.tif


Reference DEM: /Users/ben/Dropbox/UW_Shean/COP/COP30_utqiagvik_lzw-adj_proj.tif


ASP DEM: /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/stereo/20220417_2252_1040010074793300_1040010075633C00-DEM_1m.tif

Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/00.png
Plotting DEM results. This can take a minute for large inputs.
WARNING:asp_plot.stereo:

Found a DEM of difference: /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/stereo/20220417_2252_1040010074793300_1040010075633C00-DEM_1m_COP30_utqiagvik_lzw-adj_proj_diff.tif.
Using that for difference map plotting.


Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/01.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/02.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/03.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/04.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/05.png

ICESat-2 ATL06 request processing for: all
{'poly': [{'lon': -156.83342896997698, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.26658548886508}], 'res': 20, 'len': 40, 'ats': 20, 'maxi': 6, 'samples': {'esa_worldcover': {'asset': 'esa-worldcover-10meter'}}, 'cnf': 'atl03_high', 'srt': -1, 'cnt': 10}
Existing file found, reading in: atl06sr_all.parquet
Filtering ATL06-SR all

ICESat-2 ATL06 request processing for: ground
{'poly': [{'lon': -156.83342896997698, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.26658548886508}], 'res': 20, 'len': 40, 'ats': 20, 'maxi': 6, 'samples': {'esa_worldcover': {'asset': 'esa-worldcover-10meter'}}, 'cnf': 'atl03_low', 'srt': -1, 'cnt': 5, 'atl08_class': 'atl08_ground'}
Existing file found, reading in: atl06sr_ground.parquet
Filtering ATL06-SR ground

Filtering ATL06 with 15 day pad, 45 day, 91 day pad, and seasonal pad around 2022-04-17 22:52:18.495475 for: all

Filtering ATL06 with 15 day pad, 45 day, 91 day pad, and seasonal pad around 2022-04-17 22:52:18.495475 for: ground

icesat_minus_dem not found in ATL06 dataframe: all. Running differencing first.

Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/06.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/07.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/08.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/09.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/10.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/11.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/12.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/13.png


Report saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/stereo/asp_plot_report_WV03_20220417_1040010074793300_1040010075633C00_20251024_190509.pdf

%load_ext autoreload
%autoreload 2

from asp_plot.processing_parameters import ProcessingParameters

processing_parameters = ProcessingParameters(
    processing_directory=directory,
    bundle_adjust_directory=ba_directory,
    stereo_directory=stereo_directory
)
processing_parameters_dict = processing_parameters.from_log_files()

print(f"Processed on: {processing_parameters_dict['processing_timestamp']}\n")

print(f"Reference DEM: {processing_parameters_dict['reference_dem']}\n")

print(f"Bundle adjustment ({processing_parameters_dict['bundle_adjust_run_time']}):\n")
print(processing_parameters_dict["bundle_adjust"])

print(f"\nStereo ({processing_parameters_dict['stereo_run_time']}):\n")
print(processing_parameters_dict["stereo"])

print(f"\nPoint2dem ({processing_parameters_dict['point2dem_run_time']}):\n")
print(processing_parameters_dict["point2dem"])

Processed on: 2024-04-14 17:55:43

Reference DEM: /Users/ben/Dropbox/UW_Shean/COP/COP30_utqiagvik_lzw-adj_proj.tif

Bundle adjustment (0 hours and 7 minutes):

bundle_adjust -t dg --weight-image /nobackup/bpurint1/data/utqiagvik/WV/utqiagvik_wv_EE/2022/utqiagvik_10m_UTM4N_seaice_mask_0and1.tif --datum WGS84 --individually-normalize --normalize-ip-tiles --ip-per-tile 50 --matches-per-tile 10 --min-triangulation-angle 10 --mapproj-dem /Users/ben/Dropbox/UW_Shean/COP/COP30_utqiagvik_lzw-adj_proj.tif --propagate-errors --tri-weight 0.1 --tri-robust-threshold 0.1 --camera-weight 0 1040010074793300.r100.tif 1040010075633C00.r100.tif 1040010074793300.r100.xml 1040010075633C00.r100.xml -o ba/ba_50ips_10matches_dg_weight_image --threads 28

Stereo (3 hours and 41 minutes):

stereo --stereo-algorithm asp_mgm --corr-kernel 7 7 --subpixel-kernel 15 15 --cost-mode 4 --subpixel-mode 9 --corr-max-levels 5 --filter-mode 1 --erode-max-size 0 --individually-normalize --corr-memory-limit-mb 5000 --sgm-collar-size 256 --corr-tile-size 1024 --alignment-method none --corr-seed-mode 1 --compute-point-cloud-center-only --threads 24 1040010074793300_ortho_0.35m.tif 1040010075633C00_ortho_0.35m.tif ba/ba_50ips_10matches_dg_weight_image-1040010074793300.r100.adjusted_state.json ba/ba_50ips_10matches_dg_weight_image-1040010075633C00.r100.adjusted_state.json stereo_ba_50ips_10matches_dg_weight_image__ortho_0.35m_mode_asp_mgm_spm_9_corr_7_rfne_15_cost_4_refdem_COP30/20220417_2252_1040010074793300_1040010075633C00 /Users/ben/Dropbox/UW_Shean/COP/COP30_utqiagvik_lzw-adj_proj.tif

Point2dem (0 hours and 26 minutes):

point2dem --nodata-value -9999 --t_srs EPSG:32604 --threads 24 --propagate-errors --remove-outliers --remove-outliers-params 75.0 3.0 --errorimage --tr 1 -o stereo_ba_50ips_10matches_dg_weight_image__ortho_0.35m_mode_asp_mgm_spm_9_corr_7_rfne_15_cost_4_refdem_COP30/20220417_2252_1040010074793300_1040010075633C00_1m stereo_ba_50ips_10matches_dg_weight_image__ortho_0.35m_mode_asp_mgm_spm_9_corr_7_rfne_15_cost_4_refdem_COP30/20220417_2252_1040010074793300_1040010075633C00-PC.tif

import os
from asp_plot.scenes import ScenePlotter

# Optional: create directory to save plots

# plots_directory = os.path.join(directory, "asp_plots")
# os.makedirs(plots_directory, exist_ok=True)

plotter = ScenePlotter(
  directory,
  stereo_directory,
  title="Input Scenes"
)

plotter.plot_scenes(
  # Optional parameters to set for saving (requires plots_directory to be created above):

  # save_dir=plots_directory,
  # fig_fn="stereo_scenes.png"
)

from asp_plot.stereo_geometry import StereoGeometryPlotter

geometry_plotter = StereoGeometryPlotter(
  directory
)

geometry_plotter.dg_geom_plot(
    # Optional parameters to set for saving (requires plots_directory to be created above):

    # save_dir=plots_directory,
    # fig_fn="scene_geometry.png"
)

import contextily as ctx
from asp_plot.bundle_adjust import ReadBundleAdjustFiles, PlotBundleAdjustFiles

map_crs = "EPSG:32604"  # UTM Zone 4N

ctx_kwargs = {
    "crs": map_crs,
    "source": ctx.providers.Esri.WorldImagery,
    "attribution_size": 0,
    "alpha": 0.5,
}

ba_files = ReadBundleAdjustFiles(directory, ba_directory)
resid_initial_gdf, resid_final_gdf = ba_files.get_initial_final_residuals_gdfs(residuals_in_meters=True)
resid_mapprojected_gdf = ba_files.get_mapproj_residuals_gdf()
resid_triangulation_uncert_df = ba_files.get_propagated_triangulation_uncert_df()

resid_triangulation_uncert_df

plotter = PlotBundleAdjustFiles(
  [resid_initial_gdf, resid_final_gdf],
  lognorm=True,
  title="Bundle Adjust Initial and Final Residuals (Log Scale)"
)

plotter.plot_n_gdfs(
    column_name="mean_residual",
    cbar_label="Mean residual (px)",
    map_crs=map_crs,
    **ctx_kwargs
)

plotter.plot_n_gdfs(
    column_name="mean_residual_meters",
    cbar_label="Mean residual (m)",
    map_crs=map_crs,
    **ctx_kwargs
)

plotter.lognorm = False
plotter.title = "Bundle Adjust Initial and Final Residuals (Linear Scale)"

plotter.plot_n_gdfs(
    column_name="mean_residual",
    cbar_label="Mean residual (px)",
    common_clim=False,
    map_crs=map_crs,
    **ctx_kwargs
)

plotter.plot_n_gdfs(
    column_name="mean_residual_meters",
    cbar_label="Mean residual (m)",
    common_clim=False,
    map_crs=map_crs,
    **ctx_kwargs
)

plotter.title = "Bundle Adjust Initial and Final Residuals (Common Linear Scale)"

plotter.plot_n_gdfs(
    column_name="mean_residual",
    cbar_label="Mean residual (px)",
    map_crs=map_crs,
    **ctx_kwargs
)

plotter = PlotBundleAdjustFiles(
  [resid_final_gdf],
)

plotter.title = "Bundle Adjust Residual Number of Observations"

plotter.plot_n_gdfs(
    column_name="num_observations",
    cbar_label="Number of Observations",
    map_crs=map_crs,
    **ctx_kwargs
)

plotter.title = "Bundle Adjust Residuals Height Above Datum"

plotter.plot_n_gdfs(
    column_name="height_above_datum",
    cbar_label="Height above datum (m)",
    map_crs=map_crs,
    **ctx_kwargs
)

plotter = PlotBundleAdjustFiles(
  [resid_mapprojected_gdf],
  title="Bundle Adjust Midpoint distance between\nfinal interest points projected onto reference DEM",
)

plotter.plot_n_gdfs(
    column_name="mapproj_ip_dist_meters",
    cbar_label="Interest point distance (m)",
    map_crs=map_crs,
    **ctx_kwargs
)

geodiff_initial_gdf, geodiff_final_gdf = ba_files.get_initial_final_geodiff_gdfs()

plotter = PlotBundleAdjustFiles(
  [geodiff_initial_gdf, geodiff_final_gdf],
  lognorm=False
)

plotter.title = "Bundle Adjust Initial and Final Geodiff vs. Reference DEM"

plotter.plot_n_gdfs(
    column_name="height_diff_meters",
    cbar_label="Height difference (m)",
    map_crs=map_crs,
    cmap="RdBu",
    clim=(-1, 1),
    **ctx_kwargs
)

plotter.title = "Bundle Adjust Initial and Final Geodiff vs. Reference DEM (Auto Scale)"

plotter.plot_n_gdfs(
    column_name="height_diff_meters",
    cbar_label="Height difference (m)",
    map_crs=map_crs,
    cmap="RdBu",
    symm_clim=True,
    **ctx_kwargs
)

from asp_plot.stereo import StereoPlotter

# Reference DEM can be passed to the StereoPlotter, but if it is not, it will
# be searched for in the stereo logs. If it's not found there, you will be warned,
# and no difference plots will be generated.

# reference_dem = "/Users/ben/Dropbox/UW_Shean/COP/COP30_utqiagvik_lzw-adj_proj.tif"

plotter = StereoPlotter(
  directory, 
  stereo_directory,
  # Optional args:
  # dem_fn="my-custom-dem-name.tif",
  dem_gsd=1,
  # reference_dem=reference_dem,
)

Reference DEM: /Users/ben/Dropbox/UW_Shean/COP/COP30_utqiagvik_lzw-adj_proj.tif


ASP DEM: /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/stereo/20220417_2252_1040010074793300_1040010075633C00-DEM_1m.tif

plotter.title = "Hillshade with details"

plotter.plot_detailed_hillshade(
  intersection_error_percentiles=[16, 50, 84],
  subset_km=1,
)

plotter.title="Stereo Match Points"

plotter.plot_match_points()

plotter.title = "Disparity (meters)"

plotter.plot_disparity(
  unit="meters",
  quiver=True,
)

plotter.title = "Disparity (pixels)"

plotter.plot_disparity(
  unit="pixels",
  quiver=True,
)

plotter.title = "Stereo DEM Results"

plotter.plot_dem_results()

Plotting DEM results. This can take a minute for large inputs.

WARNING:asp_plot.stereo:

Found a DEM of difference: /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/stereo/20220417_2252_1040010074793300_1040010075633C00-DEM_1m_COP30_utqiagvik_lzw-adj_proj_diff.tif.
Using that for difference map plotting.

# You can also use custom color limits for the DEM results plots:

# plotter.title = "Stereo DEM Results (forced color limits)"

# plotter.plot_dem_results(
#   el_clim=(0, 10),
#   ie_clim=(0, 0.2),
#   diff_clim=(-5, 5),
# )

import glob
from asp_plot.altimetry import Altimetry

dem_fn = glob.glob(os.path.join(directory, "stereo*/*DEM_1m.tif"))[0]
aligned_dem_fn = glob.glob(os.path.join(directory, "stereo*/*DEM_1m*pc_align*.tif"))[0]

map_crs = "32604"  # UTM Zone 4N

ctx_kwargs = {
    "crs": f"EPSG:{map_crs}",
    "source": ctx.providers.Esri.WorldImagery,
    "attribution_size": 0,
    "alpha": 0.5,
}

icesat = Altimetry(
  directory=directory,
  dem_fn=dem_fn,
  # If the aligned_dem_fn is not provided, it will be calculated automatically
  # using ASP's pc_align tool, if it is available in your PATH.
  aligned_dem_fn=aligned_dem_fn
)

icesat.request_atl06sr_multi_processing(
    save_to_parquet=True,
)

# There are many more custom processing options available.
# See SlideRule's Parameter documentation for details:
# https://slideruleearth.io/web/rtd/user_guide/icesat2.html#parameters

# icesat.request_atl06sr_multi_processing(
#     res=10,
#     len=20,
#     ats=20,
#     cnt=5,
#     maxi=5,
#     save_to_parquet=True,
#     processing_levels=["ground"],
# )

ICESat-2 ATL06 request processing for: all
{'poly': [{'lon': -156.83342896997698, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.26658548886508}], 'res': 20, 'len': 40, 'ats': 20, 'maxi': 6, 'samples': {'esa_worldcover': {'asset': 'esa-worldcover-10meter'}}, 'cnf': 'atl03_high', 'srt': -1, 'cnt': 10}
Existing file found, reading in: atl06sr_all.parquet
Filtering ATL06-SR all

ICESat-2 ATL06 request processing for: ground
{'poly': [{'lon': -156.83342896997698, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.26658548886508}], 'res': 20, 'len': 40, 'ats': 20, 'maxi': 6, 'samples': {'esa_worldcover': {'asset': 'esa-worldcover-10meter'}}, 'cnf': 'atl03_low', 'srt': -1, 'cnt': 5, 'atl08_class': 'atl08_ground'}

WARNING:sliderule.sliderule:Using simplified polygon (for CMR request only!), 5 points using tolerance of 0.0001

Filtering ATL06-SR ground

ICESat-2 ATL06 request processing for: canopy
{'poly': [{'lon': -156.83342896997698, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.26658548886508}], 'res': 20, 'len': 40, 'ats': 20, 'maxi': 6, 'samples': {'esa_worldcover': {'asset': 'esa-worldcover-10meter'}}, 'cnf': 'atl03_medium', 'srt': -1, 'cnt': 5, 'atl08_class': 'atl08_canopy'}

WARNING:sliderule.sliderule:Using simplified polygon (for CMR request only!), 5 points using tolerance of 0.0001

Filtering ATL06-SR canopy

ICESat-2 ATL06 request processing for: top_of_canopy
{'poly': [{'lon': -156.83342896997698, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.26658548886508}], 'res': 20, 'len': 40, 'ats': 20, 'maxi': 6, 'samples': {'esa_worldcover': {'asset': 'esa-worldcover-10meter'}}, 'cnf': 'atl03_medium', 'srt': -1, 'cnt': 5, 'atl08_class': 'atl08_top_of_canopy'}

WARNING:sliderule.sliderule:Using simplified polygon (for CMR request only!), 5 points using tolerance of 0.0001

Filtering ATL06-SR top_of_canopy

print("There are", icesat.atl06sr_processing_levels["all"].shape[0], "ATL06SR points available.")
print("There are", icesat.atl06sr_processing_levels_filtered["all"].shape[0], "filtered ATL06SR points available.")

There are 163350 ATL06SR points available.
There are 145856 filtered ATL06SR points available.

# Filter by ESA WorldCover to remove water points
icesat.filter_esa_worldcover(filter_out="water")

print("After filtering out water, there are", icesat.atl06sr_processing_levels_filtered["ground"].shape[0], "filtered ground ATL06SR points available.")

After filtering out water, there are 45309 filtered ground ATL06SR points available.

# Apply temporal filters around acquisition date (15/45/91 day windows + seasonal)
icesat.predefined_temporal_filter_atl06sr()

# You can also create user defined temporal filters with `generic_temporal_filter_atl06sr`:

# icesat.generic_temporal_filter_atl06sr(
#     select_years=[2021, 2022, 2023],
#     select_months=[1, 2, 3, 10, 11, 12],
#     select_days=[1, 2, 3, 4, 5, 6, 7, 27, 28, 29, 30, 31]
# )

Filtering ATL06 with 15 day pad, 45 day, 91 day pad, and seasonal pad around 2022-04-17 22:52:18.495475 for: all

Filtering ATL06 with 15 day pad, 45 day, 91 day pad, and seasonal pad around 2022-04-17 22:52:18.495475 for: ground

Filtering ATL06 with 15 day pad, 45 day, 91 day pad, and seasonal pad around 2022-04-17 22:52:18.495475 for: canopy

Filtering ATL06 with 15 day pad, 45 day, 91 day pad, and seasonal pad around 2022-04-17 22:52:18.495475 for: top_of_canopy

print("There are", icesat.atl06sr_processing_levels_filtered["ground_seasonal"].shape[0], "filtered ground seasonal ATL06SR points available.")

There are 12505 filtered ground seasonal ATL06SR points available.

icesat.plot_atl06sr_time_stamps(
   key="ground",
   figsize=(15, 10),
   save_dir=None,
   fig_fn=None,
   **ctx_kwargs,
)

icesat.plot_atl06sr_time_stamps(
   key="canopy",
   figsize=(15, 10),
   save_dir=None,
   fig_fn=None,
   **ctx_kwargs,
)

icesat.plot_atl06sr(
    key="ground_seasonal",
    map_crs=map_crs,
    cmap="inferno",
    plot_dem=False,
    **ctx_kwargs
)

icesat.plot_atl06sr(
    key="ground_45_day_pad",
    map_crs=map_crs,
    cmap="inferno",
    plot_dem=False,
    plot_beams=True,
    **ctx_kwargs
)

icesat.plot_atl06sr(
    key="canopy_seasonal",
    map_crs=map_crs,
    cmap="inferno",
    plot_dem=True,
    plot_beams=True,
    **ctx_kwargs
)

icesat.mapview_plot_atl06sr_to_dem(
    key="ground_15_day_pad",
    **ctx_kwargs,
)

icesat_minus_dem not found in ATL06 dataframe: ground_15_day_pad. Running differencing first.

# Now for the pc_align aligned DEM version
icesat.mapview_plot_atl06sr_to_dem(
    key="ground_15_day_pad",
    plot_aligned=True,
    **ctx_kwargs,
)

# Show a quick histogram
icesat.histogram(
    key="ground_seasonal",
    plot_aligned=True,
)

icesat.alignment_report(
    processing_level="ground",
    minimum_points=500,
    agreement_threshold=0.25,
    write_out_aligned_dem=True,
    min_translation_threshold=0.1,
    key_for_aligned_dem="ground_15_day_pad",
)

Warning: Translation components vary by more than 25.0% across temporal filters. The translation applied to the aligned DEM may be inaccurate.

X shift range: 0.881 m (mean: 0.359 m)
Y shift range: 0.965 m (mean: 0.017 m)
Z shift range: 0.336 m (mean: -0.649 m)

Writing out: /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/stereo/20220417_2252_1040010074793300_1040010075633C00-DEM_1m_pc_align_translated.tif


Wrote out ground_15_day_pad aligned DEM to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/stereo/20220417_2252_1040010074793300_1040010075633C00-DEM_1m_pc_align_translated.tif

icesat.alignment_report_df

icesat.mapview_plot_atl06sr_to_dem(
    key="ground_15_day_pad",
    plot_aligned=True,
    **ctx_kwargs,
)

# Show a quick histogram
icesat.histogram(
    key="ground_seasonal",
    plot_aligned=True,
)

# Collect only the coincident filtereded data again for profile plotting
# icesat.filter_atl06sr(
#     h_sigma_quantile=0.95,
#     mask_worldcover_water=True,
#     save_to_csv=False,
#     select_months=[4],
#     select_years=[2022],
# )

# icesat.plot_atl06sr(
#     title=f"Cleaned beam strengths (n={icesat.atl06sr_filtered.shape[0]})",
#     filtered=True,
#     plot_beams=True,
#     plot_dem=False,
#     map_crs=map_crs,
#     **ctx_kwargs
# )

# icesat.plot_atl06sr_dem_profiles(title="Profiles", only_strong_beams=True)

from asp_plot.csm_camera import csm_camera_summary_plot

map_crs = "32619" # UTM 19S (for Uyuni)
title = "Jitter correction (Uyuni), Less Constrained"

ctx_kwargs = {
    "crs": f"EPSG:{map_crs}",
    "source": ctx.providers.Esri.WorldImagery,
    "attribution_size": 0,
    "alpha": 0.5,
}

# First set of cameras (scene 1)
original_camera = "../../tests/test_data/jitter/uyuni/csm-104001001427B900.r100.adjusted_state.json"
optimized_camera = "../../tests/test_data/jitter/uyuni/jitter_solved_run-csm-104001001427B900.r100.adjusted_state.json"
cam1_list = [original_camera, optimized_camera]

# Second set of cameras (scene 2)
original_camera = "../../tests/test_data/jitter/uyuni/csm-1040010014761800.r100.adjusted_state.json"
optimized_camera = "../../tests/test_data/jitter/uyuni/jitter_solved_run-csm-1040010014761800.r100.adjusted_state.json"
cam2_list = [original_camera, optimized_camera]

csm_camera_summary_plot(
    cam1_list,
    cam2_list,
    map_crs=map_crs,
    title=title,
    trim=False,
    figsize=(20, 15),
    add_basemap=True,
    **ctx_kwargs
)

csm_camera_summary_plot(
    cam1_list,
    cam2_list,
    map_crs=map_crs,
    title=title,
    trim=True,
    shared_scales=True,
    log_scale_positions=True,
    log_scale_angles=True,
    upper_magnitude_percentile=95,
    figsize=(20, 15),
    add_basemap=True,
    **ctx_kwargs
)

csm_camera_summary_plot(
    cam1_list,
    map_crs=map_crs,
    title=title,
    trim=True,
    shared_scales=True,
    log_scale_positions=True,
    log_scale_angles=True,
    upper_magnitude_percentile=95,
    figsize=(20, 15),
    add_basemap=True,
    **ctx_kwargs
)

	key	p16_beg	p50_beg	p84_beg	p16_end	p50_end	p84_end	x_shift	y_shift	z_shift	translation_magnitude
0	ground	0.393706	0.635813	0.840236	0.046457	0.153485	0.336891	0.223944	-0.068757	-0.593812	0.638350
1	ground_15_day_pad	0.641055	0.773392	0.921410	0.028533	0.094267	0.210883	0.155468	0.041450	-0.763709	0.780474
2	ground_45_day_pad	0.585919	0.759794	0.914996	0.033228	0.110376	0.247181	0.984863	0.071390	-0.477306	1.096756
3	ground_91_day_pad	0.560620	0.747609	0.931639	0.037673	0.127983	0.280320	0.103808	-0.461917	-0.813125	0.940912
4	ground_seasonal	0.566514	0.732346	0.886025	0.032193	0.108120	0.234966	0.324831	0.503479	-0.597760	0.846359

Comprehensive WorldView Example with `asp_plot`¶

Processing Overview¶

1. Full Report Generation¶

Individual Plots¶

2. Processing Parameters¶

3. Scene Plots¶

4. Scene Geometry Plots¶

5. Bundle Adjustment Plots¶

Bundle Adjustment Residuals (Log Scale)¶

Bundle Adjustment Residuals (Linear Scale)¶

Bundle Adjustment Observation Statistics¶

Map-Projected Residuals¶

Geodiff vs. Reference DEM¶

6. Stereo Processing Results¶

Hillshade with Details¶

Stereo Match Points¶

Disparity (Meters)¶

Disparity (Pixels)¶

DEM Results Summary¶

7. ICESat-2 Altimetry Validation¶

Request ATL06-SR Data¶

Filter ATL06-SR Data¶

ICESat-2 Temporal Distribution¶

ICESat-2 Spatial Distribution¶

DEM vs. ICESat-2 Comparison¶

PC Alignment Report¶

WIP: Profile plots¶

8. CSM Camera Model Analysis¶

CSM Camera Comparison (Basic)¶

CSM Camera Comparison (Advanced)¶

Single Camera Comparison¶

Comprehensive WorldView Example with asp_plot¶

Processing Overview¶

1. Full Report Generation¶

Individual Plots¶

2. Processing Parameters¶

3. Scene Plots¶

4. Scene Geometry Plots¶

5. Bundle Adjustment Plots¶

Bundle Adjustment Residuals (Log Scale)¶

Bundle Adjustment Residuals (Linear Scale)¶

Bundle Adjustment Observation Statistics¶

Map-Projected Residuals¶

Geodiff vs. Reference DEM¶

6. Stereo Processing Results¶

Hillshade with Details¶

Stereo Match Points¶

Disparity (Meters)¶

Disparity (Pixels)¶

DEM Results Summary¶

7. ICESat-2 Altimetry Validation¶

Request ATL06-SR Data¶

Filter ATL06-SR Data¶

ICESat-2 Temporal Distribution¶

ICESat-2 Spatial Distribution¶

DEM vs. ICESat-2 Comparison¶

PC Alignment Report¶

WIP: Profile plots¶

8. CSM Camera Model Analysis¶

CSM Camera Comparison (Basic)¶

CSM Camera Comparison (Advanced)¶

Single Camera Comparison¶

Comprehensive WorldView Example with `asp_plot`¶